Dynamics of trapped two-component Fermi gas: temperature dependence of the 
transition from collisionless to collisional regime 
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We develop a numerical method to study the dynamics of 
a two-component atomic Fermi gas trapped inside a harmonic 
potential at temperature T well below the Fermi temperature 
Tf- We examine the transition from the collisionless to the 
collisional regime down to T — 0.2 Tf and find good qual- 
itative agreement with the experiments of B. DeMarco and 
D.S. Jin [Phys. Rev. Lett. 88, 040405 (2002)]. We demon- 
strate a twofold role of temperature on the collision rate and 
on the efficiency of collisions. In particular we observe an 
hitherto unreported effect, namely that the transition to hy- 
drodynamic behavior is shifted towards lower collision rates 
as temperature decreases. 

PACS: 05.30.Fk, 71.10-w 

Recent experiments at JILA [1] on the collisions be- 
tween two oscillating spin-polarized components of a 
Fermi gas of 40 K atoms have shown that this setup is an 
important tool for investigating the dynamics of dilute 
quantum gases. These experiments have given evidence 
for a transition from collisionless (zero-sound) to hydro- 
dynamical (first-sound) behavior: the measured damping 
time r becomes very long at both low and large values 
of the estimated collision rate. The transition from the 
hydrodynamic to an intermediate regime has also been 
followed with decreasing temperature and the effect of 
Pauli blocking of the collisions from increasing occupa- 
tion of final states has been exhibited at low temperature. 

Various numerical experiments have addressed the dy- 
namics of a thermal cloud of bosons [2], even in the 
presence of a Bosc-Einstcin condensate [3] or of cold 
spin-polarized fermions interacting with a condensate via 
mean field [4]. Fermion molecular dynamics (FMD) has 
been developed [5] as a quasi-classical model for treating 
problems such as ion-atom collisions or formation of an- 
tiprotonic atoms. However, FMD is not adequate to deal 
with dilute fermionic systems, such as the neutral atomic 
gas under harmonic confinement realized in Ref . [1] , and 
an approach explicitly acknowledging the diluteness and 
other characteristics of such a system is needed. To our 
knowledge, the present work reports the first numerical 
study directed to the transport properties of ultracold 
Fermi gases. 

In this Letter we solve the Vlasov-Landau equations 
(VLE) for two-component fermionic Wigner distribu- 
tions. As in FMD [5] or in numerical studies of the dy- 
namics of thermal bosons [2-4] , the quantum-mechanical 



fluid is treated by a particle-dynamics approach. We pro- 
ceed along the path traced for a single spin-polarized 
Fermi gas [4], by duplicating it and introducing mean- 
field interactions and collisions between the two compo- 
nents. A major highlight of our numerical method is the 
strategy used to deal with collisional events. We develop 
a locally adaptive importance-sampling technique which 
allows us to handle collisional interactions faster than 
in standard Monte Carlo techniques by several orders of 
magnitude. 

Owing to this computational advance we are able to 
examine collisional effects down to T/Tp ~ 0.2, in a re- 
gion of temperature T well below the Fermi temperature 
Tp where Pauli blocking would normally grind the simu- 
lation to a halt because of numerical attrition problems 
(basically through a saturation of phase space resulting in 
vanishing efficiency of the Monte Carlo sampling). More- 
over, since we focus on the JILA setup [1] where the axial 
symmetry is maintained during the experiment, we can 
use an axially symmetric code with two-dimensional col- 
lisions, in which the angular degree of freedom is taken 
into account via an effective weight. The resulting code 
permits us to study the collisional properties of the two- 
component Fermi gas as functions of both temperature T 
and quantum collision rate T q . The numerical approach 
has full control over additional variables which are exper- 
imentally unaccessible: in particular, while in the exper- 
iments the classical collision rate can easily be estimated 
but the fully quantal collision rate T q remains unknown, 
in the numerical approach both classical and quantum 
collision rates can be counted step by step. Because of 
this control of the collision rates we are able to observe 
that, even if most collisions become forbidden classically 
and by the Pauli principle as temperature is lowered, the 
few collisions that can occur involve particles increasingly 
clustered around the Fermi level. The result is a kind of 
"catalytic effect" in which these few collisions suffice to 
drive the gas from the collisionless to the hydrodynamic 
regime. 

The model and its solution. Here we write the VLE 
for two interacting Fermi gases and describe the algo- 
rithm used in their numerical solution. Wc summarize 
some points which have been made in more detail in 
Ref. [4] and point out how the difficulties raised by Fermi 
statistics are handled in our method. 

The two fermionic components in external poten- 
tials v}H(r) are described by the distribution functions 
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f(fi (r, p, t) with j = 1 or 2. These obey the kinetic equa- 
tions 

9tfU) + m ' Vr/0) ~ Vr?70) ' Vp/0) = Cl2 t /0) ] ^ 

where the mean-held Hartree-Fock (HF) effective poten- 
tial is U^\r, t) = Vg^(r) + gn®(r, i) with j denoting 
the species different from j. Here we have set ft = 1, 
g = 2na/m r with a being the s-wave scattering length be- 
tween two atoms and m r the reduced mass, and (r, t) 
is the spatial density given by integration of f^\r,p,t) 
over momentum. 

Collisions between atoms of the same spin can be ne- 
glected at low temperature, so that in Eq. (1) the term 
C12 involves only collisions between particles with differ- 
ent polarizations: 

C 12 [f^}^ 2(2^ g 2 /V 3 
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With /« =jV)(T,p,t), f U) = 1 - fP = 

f^(r,pi,t), = 1 — f^K V is the volume occupied 
by the gas and the factors A p and A e are the usual delta 
functions accounting for conservation of momentum and 
energy, with the energies given by p 2 /2rrij + U^(r,t). 

The numerical procedure by which the VLE is ad- 
vanced in time consists of three basic steps: (i) initial- 
ization of the fcrmionic distributions, (ii) propagation in 
phase space, and (in) collisional interactions. 

The initial distributions in equilibrium at the bottom 
of a harmonic trap are generated by using the HF ex- 
pression 
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where (3 — l/fc^T and ^ is the chemical potential 
of species j [6]. The particle densities entering U^\r) 
are to be determined self-consistently by integration over 
momenta, and the momentum distributions of the two 
clouds need generating. 

To exploit the axial symmetry of the system, we 
hrst dchnc the angularly integrated particle densities 
AfW> (r, z) = 2nrn^ (r, z) on a {r, z} grid in cylindrical 
coordinates and we move to a particles-in-cell descrip- 
tion by locating a number Af^ (r, z)ArAz of fermions 
in each cell of volume ArAz [4]. A momentum distri- 
bution with low statistical noise is generated by rep- 
resenting each fermion by means of N q computational 
particles ("quarks"). The i-th quark is located at point 
{Pir > P%e ) Viz } hi momentum space by using a Monte Carlo 
sampling and making sure that each cell of volume h 3 is 
occupied by no more than N q quarks. This control in 3D 
phase space is transferred to 2D by imposing a maximum 



number wN q of quarks in the 2D cell {Ar, Az, Ap ri Ap z } 
of volume h 2 . Here the weight w = int(2pi?r/ft) takes into 
account that the number of available cells in 2D depends 
on the radial position and on the number of particles 
through the Fermi momentum pp = ^/2mhkBTp. 

In the propagation step the two clouds are rigidly dis- 
placed from the center of the trap along the z direction 
and start evolving in time in the {r, z) plane by per- 
forming oscillations at their respective frequencies. We 
exploit the fact that by symmetry the average value of pe 
for each value of r is zero and does not change in time, 
to decouple the variable pg from the equations of motion. 
We make the approximation that the angular momentum 
of each quark is left unchanged by the collisional events 
and take into account the third dimension in the exclu- 
sion principle by suppressing collisions whose hnal states 
are occupied by more than wN q quarks. At each time 
step the quarks are moved by the confinement and the 
mean-held forces by using a Verlet algorithm on a grid 
of mesh spacing dx > v dt, where v is a typical particle 
velocity and dt the time step. This inequality is dictated 
by accuracy and stability of the propagation step [4] . No 
exclusion constraint is applied during this step: this does 
not lead to any appreciable deviation from Fermi statis- 
tics down to 0.2 TV, where the system is still sufficiently 
dilute. This has been checked by monitoring violations 
of the Pauli principle at each time step. 

Finally we come to the collision step, which involves 
most of the innovative aspects of our scheme. Collisions 
are tracked on a grid of mesh spacing 5x of the order 
of the de Broglie wavelength As, which is smaller than 
the particle mean free path I and larger than the spacing 
dx of the propagation mesh (I > Sx ~ \b > dx). The 
hrst inequality enhances the statistical accuracy of the 
collision step, whereas coarse graining with respect to 
the propagation step avoids the need for Pauli-principlc 
constraints, at least down to 0.2 Tp. 

Before turning to the results we add some technical 
details on how a speed-up of the code by several orders 
of magnitude has been achieved. The number of prob- 
able collisions between all possible pairs of quarks be- 
longing to the two species in each cell of volume h 2 is 
evaluated at each computational time step as dN co u = 
dt\J Ni N2 v ij Cij ■ Here Vij is the magnitude of the 
relative speed of particle i of species 1 and particle j of 
species 2, cry is the corresponding differential cross sec- 
tion, and N\ (N%) is the number of particles of species 1 
(2) in the given spatial cell. If dN co u < 1, the collision 
probability is accumulated over the subsequent time steps 
until an integer number dI co u = int(d/V co ;z) of collisions 
occurs. Within each spatial cell the pairs of particles are 
ordered according to the value of their classical collision 
probability, from largest to smallest. The acceptance rate 
of the MonteCarlo sampling is enhanced by two-three or- 
ders of magnitude at each step by hltering out pairs with 
classical probability smaller than a given threshold. This 
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is a form of importance sampling and the threshold is dy- 
namically adjusted cell-by-cell in such a way as to guar- 
antee the correct supply of dI co u collisions at each time 
step. The collision probability becomes smaller than the 
classical one after multiplication by the quantum sup- 
pression factor 1 — N(r, z,p r ,p z )/(wN q ) due to the occu- 
pancy of the final state, and consequently the selection of 
the most likely pairs to collide becomes essential. To this 
purpose each particle is allowed to collide only with the 
partner which maximizes the product VijUij. The origi- 
nal pool of N\N 2 collisions is cut down to N\ collisions 
only, with a resulting additional speed-up of at least one 
order of magnitude. 

Results. As a first application of the numerical 
method we have considered 200 magnetically trapped 
40 K atoms, which are represented by a total number of 
4x 10 3 quarks. As in the JILA experiments [1], the atoms 
are equally shared among two different spin states (m/ = 
9/2 and mj = 7/2) in harmonic traps with slightly dif- 
ferent longitudinal frequencies (wg/ 2 = 2-7rxl9.8 s _1 and 
^7/2 = 27TX17.46 s _1 ). When the two clouds after initial- 
ization are rigidly displaced from the center of the traps, 
in the absence of interactions they would keep oscillating 
at their respective trap frequencies without damping. 

The collision rate T q is independently varied by chang- 
ing the scattering length a, thus mimicking the exploita- 
tion of a suitable Feshbach resonance [7]. To offset 
the difference between the number of atoms used in 
the simulation (TV = 200) and that in the experiments 
(N exp <~ 10 6 ), the reference value of the off-resonant scat- 
tering length is scaled by a factor (Nexp/N) 1 / 2 <~ 10 2 , 
thus producing a system with essentially the same col- 
lision rate and mean field potential as in Ref [1]. The 
transition from the collisionless to the collisional regime, 
as driven by varying the scattering length at various tem- 
peratures, is shown both in the plot of the frequency of 
the dipole mode for the two components in Fig. 1 (a) and 
in the plot of the damping rate 7 = 1/t of the axial 
motion of the centers of mass Zcm(i) in Fig. 2. 

In Fig. 1(a) the oscillation frequencies utj at T = Tp 
and T = 0.2 Tp, with k B T F = ftw 9/2 (67V 9/2 ) 1 / 3 , are eval- 
uated by fitting ZcH(t) with the functions z cos(cjji)e~ 7t . 
At very low T q the dipole mode frequencies are given by 
the corresponding trap frequencies with a shift due to 
the mean-field potential. At intermediate values of T q 
the data points show large fluctuations, due to the fact 
that in this region just a few collisions can drastically al- 
ter the motion of the clouds. However, the trend towards 
a locking of the two dipole mode frequencies at large T q 
is very clear and the location of the locking is identified 
with reasonable accuracy. 

The main new physical result of this study is the shift 
of the locking transition to lower T q as temperature is de- 
creased. This is shown in Fig. 1(b). This effect is related 
to Fermi statistics: at low temperature the collisions in- 



volve particles on a narrower region around the Fermi 
level and have a greater impact on the global dynamics 
of the gas. A smaller number of collisions is needed to 
produce locking of the two interacting species. 

The damping rate 7 obtained from the fit of Zcm(t) is 
essentially the same for the two components and is shown 
in Fig. 2. We have also evaluated the correlation function 
(j){t) =< \Z^(t')\\Z^ 3 2(t' + t)\ > between the magnitude 
of the turning points i?cm, which decays exponentially 
as exp (—7*) at large t. This estimate of 7 is in good 
agreement with that obtained from the center-of-mass 
motions, but yields a less noisy signal in the intermedi- 
ate region. In the collisionless regime the damping rate 
increases linearly with T qi while in the collisional one it 
scales like T^ 1 . This is seen in Fig. 2, which also shows 
again that the hydrodynamic regime is reached at lower 
T q as T is lowered. 

The collision rate T q can be driven by cooling at fixed 
scattering length, as is seen from Fig. 3. The various 
curves, after scaling by a factor proportional to the clas- 
sical cross-section, show a residual weak dependence on 
a, which is due to the mean-field interaction between 
the two overlapping fermionic components. The trend of 
these curves is mostly classical and only at very low tem- 
perature (T =0.3-0.2 Tp) Pauli blocking becomes mani- 
fest, as signalled by the collapse of the curves into a single 
one. Our results thus suggest new experiments, in which 
the collision rate is changed via a thermal drive. 

The implication of the results shown in Fig. 2 and 
Fig. 3 is that this thermal drive rests upon the increas- 
ing importance of a decreasing set of "strategic" collisions 
involving particles in states clustered around the Fermi 
level. Our data also indicate that thermal cooling will 
need to reach temperatures below 0.2 Tp in order to see 
the complete transition from the collisional to the col- 
lisionless regime. Efforts to develop a concurrent code, 
including the Pauli principle in the Lagrangian evolution 
to treat the collisional Fermi gas well below 0.2 Tp, are 
currently under way. 
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FIG. 3. The collision rate r 9 (in units of s 



scaled 



by the factor 10 8 (ao/a) 2 with ao being the Bohr radius, 
as a function of temperature T (in units of TV) for vari- 
ous values of the scattering length a. From bottom to top: 
a = (12, 15, 18, 21, 24, 27, 30, 33) x 10 3 a . 




(b) 

FIG. 1. (a) The oscillation frequencies lo (in units of s _1 ) as 
functions of the quantum collision rate F q (in units of s _1 ) for 
the two components of the gas at T = TV ( x ) and T — 0.2 Tf 
(+). The horizontal dashed lines show the bare trap frequen- 
cies, (b) The collision rate F q at frequency locking as a func- 
tion of temperature T (in units of TV). 
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FIG. 2. The damping coefficient 7 (in units of s _1 ) as a 
function of the collision rate F q (in units of s _1 ) for T — Tf 
(filled squares), T = 0.67V (crosses) and T = 0.27V (empty 
squares). In the collisionless region 7 has been fitted by the 
function j(T q ) = aiT q and in the collisional regime by the 
function 7(r g ) = c^/r,: the fits are shown by a continuous 
line for T = Tf, a dashed line for T = 0.67V and a dotted 
line for T = 0.27V. 
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